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ABSTRACT 

With the ever increasing resolution of A-body simulations, accurate subhalo detection 
is becoming essential in the study of the formation of structure, the production of 
merger trees and the seeding of semi-analytic models. To investigate the state of halo 
finders, we compare two different approaches to detecting subhaloes; the first based on 
overdensities in a halo and the second being adaptive mesh refinement. A set of stable 
mock NFW dark matter haloes were produced and a subhalo was placed at different 
radii within a larger halo. SUBFIND (a Friends-of-Friends based finder) and AHF (an 
adaptive mesh based finder) were employed to recover the subhalo. As expected, we 
found that the mass of the subhalo recovered by SUBFIND has a strong dependence 
on the radial position and that neither halo finder can accurately recover the subhalo 
when it is very near the centre of the halo. This radial dependence is shown to be 
related to the subhalo being truncated by the background density of the halo and 
originates due to the subhalo being defined as an ovcrdcnsity. If the subhalo size is 
instead determined using the peak of the circular velocity profile, a much more stable 
value is recovered. The downside to this is that the maximum circular velocity is a poor 
measure of stripping and is affected by resolution. For future halo finders to recover 
all the particles in a subhalo, a search of phase space will need to be introduced. 

Key words: methods: numerical - galaxies: formation - galaxies: haloes - cosmology: 
theory - dark matter 



1 INTRODUCTION 

It has long been understood t hat dark matter plays an essen- 
tial role in galaxy formation. IWhite £T Rocs ( 197i|) demon- 
strated that dark matter haloes act as potential wells within 
which infalling material can be captured and condense to 
form galaxies. As the universe ages, these haloes merge to 
form larger structures and this continued process produces 
the framework of the universe that we see today. This so 
called hierarchical model of galaxy formation has been put 
to many tests including those generated by A-body simu- 
lation. One of the most wid ely used of these sim ulations is 
the Millennium Simulation (|Springel et al.| [2005) which ac- 
curately reproduced the large-scale structure of a 500 Mpc/h 
cube region of the universe. 

One of the challenges of studying the results of A-body 
simulations has been finding a consistent way of identifying 
the structures and substructures within them. Detailed stud- 
ies of haloes and subhaloes require halo finders, codes that 
scan the simulation outputs and identify structures. Many 
different halo finders are available and each uses different 
techniques and definitions of the haloes they find. Broadly, 
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halo finders fall into two general categories; those based on 
the Friends-of-Friends (FoF) technique and those based on 
grids. 

FoF was first proposed by iDavis et al.l (|l985t ) and lo- 
cates haloes based on a predetermined linking length for 
particles. This is usually a fraction of the mean inter- 
particle separation and any two particles closer than this 
distance are linked together. Isolated sets of linked parti- 
cles are then identified as the haloes. Commonly a value 
of 0.2 times the mean inter-particle separatio n is chosen 
motiv ated by SCDM (fi = 10 & fi A = 0.0) (jPavis et al.1 
Il985h and a slightly lower value of O.f 6 i s sometimes adopted 
for ACDM (Q = 0.3 fc 0\ = 0.7) |Lacev fc Cole! 1 19931 : 
lEke. Cole fc Frenkl fl996h . Despite the difference, conver- 
gence between cosmologi es in the halo mass function can 
be found using 0.2 (see I Jenkins et al.l 1200 ll ) making this 
the most widely used. The FoF method was im plemented 
in, for example, su bfind l|Springel et al.l 120011 ) and hfof 
l|Klypin et al.1 If 9991) , with different techniques being used 
to find subhaloes. hfof uses hierarchical FoF to locate the 
subhaloes by using a shorter linking length inside the halo, 
while subfind searches the haloes for overdensities in the 
density profile. 

The grids method of halo finding works by placing a 
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grid across the simulation and smoothing the discrete par- 
ticle data onto that grid and then locating the densest cells. 
Refinement can be built onto the grid to obtain improved 
resolution and to increase the speed of the code. The den- 
sity peaks that are located on the grid can then be used 
as the seeds for potential structures. This technique was 
used by, for example, ahf i Knollmann fc Kne"bell2009l ) and 
ASOHF (|Planelles fc Quilia 120101 1. The variations between 
these codes comes in the definition of haloes. AHF uses iso- 
density contours on the grid, while ASOHF uses spherical 
overdensities. 

FoF and grid based methods are the two main 
ways for locating structure, bu t there are alternatives . 
More recent finders, such as hsf (|Macieiewski et all [2009), 
have tried using phase space to identify subhaloes. This 
extends the search based on position and density to 
incorporate the velocity of the particles. Bulk veloci- 
ties can then also be used to help identify structures. 
Other finders that have tried different techn iques in- 
clude voboz |Nevrinck, Gnedin fc Hamilton! [20051 ). which 
replaced the uniform grid with a Voronoi diagram, and SUR V 
jTormen. Moscardini fc Yoshidal l20oi ; iGiocoli et aIll201Cf) . 
which uses knowledge of the structures from one snapshot 
to help find structure in the next. While this summary of 
halo finders is by no means exhaustive, it does give a flavour 
for the different techniques employed. A thorough review of 
the different types of halo finders available and their effec- 
tiveness will be found in Knebe et al. (in preparation). 

The importance of accurate subhalo detection has in- 
creased in recent years with the advances in high res- 
olution simulations. Various simulations of Milky Way 
sized haloes have been produced including via Lacte a 
jDiemand. Kuhlen fc Madaul 120071 ; iDiemand et~aH l2008h 
Aqua rius 1 Springel et al.l 120081 ) and ghalo ( Stadel et all 
2009). As ex pected, these hal oes contain a wealth of sub- 
structure fsee lGao et al.ll2004 ). However, it is important to 
ask how robust the recovered properties of subhaloes are 
to the choice of subhalo finder. For example, subhaloes are 
identified initially as overdensities in their host haloes. We 
expect picking out such overdensities to be more difficult 
in the innermost parts of the host haloes where the back- 
ground density is the greatest. If one halo finder is less able 
to pick out these overdensities than another halo finder, we 
would expect this halo finder to systematically underpredict 
the numbers of subhaloes in the inner parts of haloes, which 
would have important implications for how we interpret the 
results of, for example, the radial distribution of subhaloes 
and subhalo mass loss. 

In this paper we set out to quantify the extent to which 
our choice of halo finder impacts on the radial distribution of 
subhaloes that we re cover. Spec ifically we focus on subfind 
(|Springel et al.ll200ll ) and ahf (|Knollmann fc Knebdl2009T ) 
and ask how well these hal o finders can recover the prop - 
erties of a NFW subhalo (|Navarro. Frenk fc White] Il997l ) 
embedded in a more massive host NFW halo. The advan- 
tage of this approach is that, unlike using haloes and sub- 
haloes drawn from cosmological simulations, we know ex- 
actly which particles belong to the host and to the subhalo 
at initial time and we can track their positions and veloci- 
ties at all subsequent times. This provides a clean test of the 
halo finders because any discrepancies found can be identi- 
fied easily. 



The rest of this paper is setout as follows. In Sj2] we 
outline the methods used, including summaries of the halo 
finders and the process of constructing a mock 6D (x, y, z, 
Vx, v y , v z ) NFW halo by reproducing the density and veloc- 
ity profiles. We then use this construction, in |J3] to model 
an infalling subhalo. This is undertaken in two ways, first 
by considering how well the halo finders recover the subhalo 
when simply placed at different radii within the main halo. 
The second method is to let the subhalo fall into the main 
halo under gravity and compare how the different halo find- 
ers recover the subhalo. Having established the accuracy of 
the halo finders, in §Z\ we investigate the effect the trajec- 
tory of the subhalo has on stripping as it passes through the 
halo. In iJH we test the reliability of recovering the peak in 
the circular velocity profile. Finally we summarise our re- 
sults. Throughout this work, a standard ACDM cosmology 
has been adopted, taking Qo = 0.3, Q,a = 0.7 and h = 0.73, 
where appropriat e, consistent with ob servations from wmap 
first year results (|Spergel et al.ll2003l ). 



2 METHODS 
2.1 Halo Finders 

For the purpose of this work we focus on two halo finders 
that rely on different methods to detect haloes and sub- 
haloes. 

2.1.1 AHF 



Alli{3( Knollmann fc Kncbc 



2009) is an updated version of 



mhf IjGill. Knebe fc Gibsonll2004r ) and works using an adap 



tive mesh refinement method. It begins by placing a user- 
defined grid across the box and calculates the particle den- 
sity in each cell. If this is greater than a user-specified value, 
then the cell is refined with a smaller grid. The particle den- 
sity is then recalculated on this finer grid and, if required, 
further refinement is carried out. Once all the refinements 
are carried out, a hierarchical grid tree of the density dis- 
tribution has been produced and this can be used to find 
structure. Throughout this work, we used a grid of 128 cells 
with refinement being carried out in cells that contain more 
than 3 particles. 

The most refined and isolated cells are used as poten- 
tial halo centres and these are linked to the coarser grids 
to build the structure. If two isolated centres join up on 
a coarser grid then these are combined into one structure. 
By considering these separate, isolated points in one struc- 
ture, substructure can be defined. Once the structures are 
identified, starting on the lowest level of substructure, they 
are tested for boundness in isolation. This is conducted by 
comparing the particles velocity to the local escape velocity 
obtained using a spherical potential approximation. If a par- 
ticle is found to be unbound it is assigned to the next highest 
level of structure until it is dispensed with if not bound to 
the halo. The haloes are then truncated at the virial ra- 
dius (see H2.2p to define their size. For the subhaloes, not 
all have a low enough overdensity to satisfy the virial radius 
due to the background density of the halo. If this is the case 

1 Available from http://popia.ft.uam.es /AMIGA 



© 2010 RAS, MNRAS 000,[T}{9] 



The Accuracy of Subhalo Detection 3 



then they are truncated by a sharp spherical boundary at 
the outer radius at which their density profile first shows an 
upturn and starts to rise with increasing distance. 



For fi = 0.3, Q A = 0.7 and z = 0.0, A vir w 101. Using the 
scale radius and the virial approximation, the characteristic 
density is given by, 



2.1.2 SUBFIND 



5c = 



SUBFIND l|Springel et al.ll200~lh begins by conducting a stan- 
dard Friends-of-Friends (FoF) search of the simulation vol- 
ume to identify haloes. At each particle the local density 
is then calculated using a local SPH-like smoothing kernel 
interpolation over the nearest neighbours. Any locally over- 
dense region is then considered as a subhalo candidate with 
its shape being defined by an isodensity contour that tra- 
verses the saddle point in the density profile of the halo. 
This is found by lowering the global density threshold and 
selecting out the overdense regions. At this stage particles 
can be members of more than one structure allowing differ- 
ent levels of substructure to be determined. For this work, 
we used a FoF linking length of 0.2 and 10 particles for 
the SPH density calculation allowing SUBFIND to recover all 
subhaloes with 10 or more particles. Tests were also carried 
out using higher values for the SPH density calculation, but 
the number of particles recovered was found to be relatively 
insensitive to this parameter for the size of the subhalo we 
used. 

Once subhalo candidates have been identified, an un- 
binding procedure is used to determine iteratively which 
particles are not gravitationally bound. This is achieved by 
defining the centre of the subhalo as the position of the most 
bound particle and the bulk velocity as the mean velocity of 
the particles in the group. The kinetic and potential energies 
of the particles are then compared and unbound particles are 
removed. The final step is to assign particles that are listed 
in multiple structures to just one. To solve this, the parti- 
cles are assigned to the smallest structure they are found in. 
The remaining FoF particles that have not been assigned to 
substructure are then tested for boundness and assigned to 
the background halo. Any particles that are not bound to 
anything are then classified as FoF 'fuzz'. 

2.2 Constructing a Mock Halo 

The following outlines the process of constructing a mock 
dark matter halo. For simplicity we have limited ourselves 
to the case of a spherical halo that follows a radial NFW 
density profile, 

Pcrit^c 



p(r) = 



(1) 



r/r s (l + r/r s ) 2 ' 

where p cr it is the critical density of the universe, r s is the 
scale radius and S c is the characteristic density. Dark matter 
haloes are characterised by their virial mass, 



\* 47r 3 A 

= — »\dr Avir Pcrit, 



(2) 



where r v i r is the virial radius and A v j r is the virial approxi- 
mation given bv lBrvan fc Normanl ()l998l ) as, 



Avir = 18tt + 82(fi(z) - 1) - 39(H(«) - l) 2 , 



where, 



O(z) 



n (l + z) 3 
fio(l + zf + fi A ' 



(3) 



(4) 



(5) 



3 ln(l + c) - c/(l + c) ' 

where c = r v i r /r s is the concentration. 

Using these conditions, a Monte Carlo realisation can be 
constructed by defining the number of particles within r v i r , 
Avir, and specifying the concentration of the halo required. 
The Monte Carlo realisation is produced by drawing a ran- 
dom enclosed mass and inverting to find a radius. This is 
then turned into a set of coordinates by specifying they pro- 
duce a smooth distribution on the surface of a sphere. The 
mass of a NFW halo continues to increase with increasing 
radius and so in principle has infinite mass; we circumvent 
this by truncating the halo beyond a cut-off radius, r cu t- 
This modifies the density profile so that p(r < r cut ) follows 
the NFW profile and p(r > r cut ) = 0. For this work we set 
r C ut = 2r v i r . A smoother truncation could be produced by 
using a exponential decay at the edge of the halo. 

Once the halo is constructed, each particle needs to be 
given a velocity that reproduces the velocity dispersion, <r(r) , 
of a halo. Dark matter haloes are supported by the random 
motion of the particles and to get an accurate representation 
we need to reproduce this in the velocity of the particles. The 
velocity dispersion can be obtained by considering the Jeans 
equation, 



IA (pcT?) + 2 ^ 
p dr r 



dr 



(6) 



where j3 — 1 — ag(r)/a^(r) and $ is the gravitational po- 
tential. Assuming isotropy, o~g(r) = cr x {r), (3 = and the 
velocity dispersion is given by, 



°V 2 (r) = ~t~t 



p(r ) — dr . 



(7) 



P( r ) J r dr 

This integral was solved by lLokas fc Mamonl (|200ll ). and 
confirmed here, to give, 



c 2 s(l + cs) 



V 2 



(1 + cs) 2 



2[ln(l + c) - c/(l + c) 



1 + CS » n*a* 



ln(cs) 



cs 1 + cs 



x ln(l + cs) +3 In (1 + cs) + 6Li 2 (-cs) 



(8) 



where s = r/rvir, VW is the circular velocity at the virial 
radius and Li 2 (x) is the dilogarithm function given by, 



Li 2 (x) 



ln(l - t) 



dt. 



(9) 



The 3D velocity dispersion is then given by the sum of the in- 
dividual components. Since isotropy was assumed this gives 
°"3d( 7 ") — 3a 2 (r). To generate a velocity distribution func- 
tion for a given rad ius, a Maxwell- Boltzmann distribution 
can be assumed fcf. lHernquist|[l993l ). 



F(v,r) = 4tt 



( 1 \ 3 / 2 



exp 



2<j2 



(10) 



Note that the dilogarith m approximation given in equation (17) 
of lLokas fc Mamorj ifeoOlll is not suitable for this task. 
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Figure 1. The density profile of a M vir = 1O 14 M , N vir = 10 6 
and c = 5 halo left to evolve over 8 Gyr. The black line denotes 
the theoretical NFW profile, while the mock halo is shown ini- 
tially (black pluses), after 4 Gyr (red asterisks) and 8 Gyr (blue 
crosses). The arrow represents the Plummer equivalent softening 
(h = 2.8e = 8.4kpc). 



Figure 2. The fraction of particles recovered at a given separa- 
tion as the subhalo is placed at different positions within the halo. 
Both halo finders recover consistent sizes across the multiple real- 
isations, resulting in small error bars. The dotted line represents 
the fraction of particles recovered if the subhalo is truncated at 
the radius where its density is equal to the background density of 
the halo. 



The function F(v, r) is normalised such that, 




The velocity of each particle can then be obtained using the 
probability distribution of equation (|10[) . Having obtained 
the density and velocity profiles of the halo, the only thing 
left is to assign a direction to each velocity. This is done by 
simply requiring that the directional velocity vectors pro- 
duce a smooth distribution on the surface of a unit sphere. 

To test the stability of this setup, an isolated halo 
with Af vir = 1O 14 M , iV v i r = 10 6 a nd c = 5 was left to 
evolve over 8 Gyr using GADGET-2 (jSpringell [2005'). The 
spline gravitational softening was set to e = 3kpc corre- 
sponding roughly to the radius of the 100th particle (see 
IPower et al.l 120051 ). Fig. [I] shows that the halo retains the 
overall shape of an NFW profile, except at the centre 
where the profile has flattened similar t o that observed by 
iKazantzidis, Magorrian fc Moorel (|2004l ). This flattening of 
the density profile is caused by approximating the distribu- 
tio n function with a Maxwe ll-Boltzmann. As demonstrated 
in IKazantzidis et aTl |2004l '). this will lead to an over es- 
timate of any stripping that occurs. Despite this, it will 
have no effect on the ability of halo finders to recover the 
haloes. Thi s was c onfirmed by using the m ethod outl i ned in 
I Read et al. (20 0tj) to generate haloes with iPlummerl (|l91ll ) 
and lHernquistl (| 1990T ) density profiles based on their 6D dis- 
tribution functions. When the same tests were carried out 
on these haloes, the same patterns between the halo find- 
ers was found as for the NFW with the Maxwell-Boltzmann 
approximation. 



3 MODELLING AN INFALLING SUBHALO 
3.1 Static Infall 

The first method of modelling the infall of a subhalo we 
adopted was to consider how well different halo finders re- 
covered the subhalo at a given radius. This was achieved 
by placing the same sized subhalo by hand at different radii 
within the main halo and attempting to recover it with each 
halo finder. A halo was generated with M v i r = 1O 14 M0, 
7V vir = 10 B and c = 5 and a subhalo with M vil = 10 12 M Q , 
A^vir = 10 4 and c = 12. The concentration of the sub- 
halo was set to be higher than the halo in order to re- 
flect the conditions found in cosmological simulation s (see 
iBullock et all l200ll ; lEke, Navarro &; Steinmetz|[200ll ). The 
subhalo was then placed at different distances away from 
the centre of the halo and given a velocity, 

2GM halo , . 

v = \— ' ( 12 ) 

y I scp 

where Mhaio is the mass of the halo and r aep is the separation 
of the centres of the halo and subhalo, towards the centre 
of the halo. This velocity corresponds to the conversion of 
potential energy to kinetic, for two point masses, as the sub- 
halo falls in from infinity. When the subhalo was placed at 
the centre of the halo, r scp = 0.0 so v — 5* 00. To overcome 
this, the subhalo was given a velocity of the previous closest 
separation when it was at the centre of the halo. This set-up 
was produced 100 times for each separation using different 
random number seeds. Consistent realisations were found 
each time. 

Fig. [2] shows the fraction of particles recovered by each 
halo finder at different separations. Neither halo finder can 
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Figure 3. The maximum circular velocity of the recovered sub- 
halo as it is placed at different separations. Both halo finders 
accurately recover the peak, with a small radial dependence dis- 
played in SUBFIND. 



recover the subhalo when it is near the centre of the halo. 
This corresponds to the densest region of the halo and leads 
to any overdensity from the subhalo being hidden. As the 
separation is increased AHF has a steep rise in the fraction 
of particles it recovers until it is finding the complete sub- 
halo from ~ 0.5r v ir outwards. SUBFIND does not have such 
a drastic change and continues to underestimate the size of 
subhalo all the way out to ~ 1.5r v i r . 

We can gain some insight into the strong radial depen- 
dence in recovered particle number in SUBFIND by consider- 
ing the following simple argument. SUBFIND identifies sub- 
haloes as overdensities; it identifies when a subhalo's local 
density equals its host halo's local density. This equates to, 

' X ! * (13) 

1+ " 



where <5c halo and 6c Bub are the characteristic densities of the 
halo and subhalo respectively (equation [SJ , r Shalo and r Saub 
are the scale radii of the halo and subhalo respectively, r scp 
is the separation of the centres of the halo and subhalo and r 
is the radius of the subhalo at which the densities are equal. 
The number of particles within r cannot exceed iVvir by con- 
struction. The shape of the theoretical curve (dotted line in 
Fig. [2j implied by equation (|13|) reasonably captures the 
shape of the curve recovered by SUBFIND. The agreement is 
not perfect, equation (|13[1 predicts more mass should be re- 
covered at larger radii than is recovered in practise, but the 
differences can be easily understood. First, based on the ran- 
dom nature of the velocity assignment some of the particles 
will have large velocities and will therefore not be bound. 
The effect of this will be to cause the two curves to devi- 
ate systematically from each other with increasing radius. 
Second, SUBFIND identifies overdensities as saddle points in 
the mass density profile rather than by equating subhalo 



Figure 4. The fraction of particles recovered at a given radius 
as the subhalo is allowed to fall into a halo from infinity. The 
subhalo experiences the most stripping when it passes through 
the centre of the halo. Neither halo finder can detect the subhalo 
as it passes through the centre of the halo and they yield different 
sizes for the subhalo either side of this region. 



and halo mass profiles, as implied by equation (113[) . Overall 
the curve shares the same shape as that found using sub- 
find, indicating that the background density is affecting the 
ability to recover the subhalo. 

Implanting a NFW subhalo in a larger halo, defining 
the virial radius using equation (O, is obviously a highly 
idealised situation. Realistically the subhalo would be ex- 
pected to undergo stripping which would cause it to be 
stripped down to its tidal radius at different points within 
the halo. This tidal radius would roughly correspond to the 
radius at which there is a saddle p oint in the density pro- 
file l|Tormen. Diaferio fc Svei1ll998l ). This also corresponds 
to the size of the overdensity that SUBFIND is recovering. 
Therefore, if the edge of the subhalo is defined as the tidal 
radius, SUBFIND would give consistent recovery of the sub- 
halo. 

A different method of determining the size of the sub- 
halo is_to_conader_the peak in the circular velocity profile 
fsee lGhigna et alJll998ll200Ch . This will be less affected by 
truncation of the subhalo, as the particle with the maxi- 
mum circular velocity is closer to the centre. Fig. [3] shows 
the recovered maximum circular velocity for the subhalo at 
different separations. This was obtained by calculating the 
circular velocity for each particle in the subhalo and tak- 
ing the largest of these as the peak. As expected, both halo 
finders more accurately recover the subhalo size using this 
method. SUBFIND still displays a slight radial dependence, 
with a gradual decrease towards the centre of the halo. This 
is caused by high velocity particles near the centre of the 
subhalo being unbound due to the truncation. As the sub- 
halo was not detected at the centre of the halo, it is not 
possible to obtain a circular velocity there. 
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Figure 5. Fraction of particles recovered (upper panels) and maximum circular velocity (lower panels) for the subhalo as a function of 
time as the subhalo falls through the halo. For each case the subhalo is given a velocity along the x-axis toward the halo and starts offset 
by 3.0r v i r in the x-axis and 0.0 (black line), 0.5r v i r (red line) and 1.0r v i r (blue line) in the y-axis. This corresponds to a closest radial 
approach to the centre of the halo of 0.0, 0.2r v i r and 0.5r v i r respectively. 



3.2 Dynamic Infall 

The second method of investigating the infall of a subhalo 
was to allow the system to evolve under gravity. The same 
halo and subhalo properties were set up as in Q3.1\ The sub- 
halo was then placed so that r scp = 3rvi r of the halo and 
it was given a velocity toward the centre of the halo from 
equation (|12|l . The subhalo was then left to free-fall through 
the halo for 6 Gyr using GADGET-2 with gravitational soften- 
ing e = 3kpc. Snapshots were taken every 0.05 Gyr. During 
this run cosmological expansion was turned off so the haloes 
were only affected by gravity. 

Fig. [4] shows the fraction of particles recovered by SUB- 
find and AHF as the subhalo passed through the halo. 
The subhalo undergoes a large amount of stripping, loosing 
around 75 per cent of its mass. Most of this stripping occurs 
as the subhalo passes through the very centre of the halo. 
This corresponds to the greatest rate of change of the poten- 
tial and so would be expected to have the largest effect. As 
predicted in H3.1l both halo finders fail to recover the subhalo 
as it passes through the centre of the halo and disagree about 
the size of the subhalo immediately either side of this region. 
The largest discrepancy occurs when the subhalo is within 
the virial radius of the halo. As expected due to its defini- 
tion of a subhalo, SUBFIND recovers a smaller subhalo during 
the infall phase compared with AHF. After the subhalo has 
passed the centre of halo, AHF recovers a much larger num- 
ber of particles due to its unbinding procedure being less 



efficient and this is discussed further in Sj4] A s expected, the 
level o f stripping obser ved is consistent w i thlHavashi et all 
l|2003h and higher than iKazantzidis et af] (|2004h . 



4 SUBHALO STRIPPING 

As seen in H3.2I an infalling subhalo only undergoes stripping 
as it passes through the very centre of the halo. This should 
mean that any subhalo that does not pass through the centre 
of the halo and is merely deflected around it should undergo 
significantly less stripping. To test this hypothesis, the sub- 
halo was placed at a separation of 3.0r v i r in the a:-axis and 
0.0, 0.5r v ir and 1.0r v i r in the y-axis. In each case the sub- 
halo was given the same velocity along the x-axis toward 
the halo as in £13.21 The subhalo that was on the x-axis fol- 
lowed the same path as the subhalo in >i3. 21 passing straight 
through the halo centre. The other two subhaloes were de- 
flected around the halo centre with closest approaches of 
0.2r\,ir and 0.5r v i r respectively. 

Fig. [5] shows the fraction of particles recovered by each 
halo finder for the three scenarios outlined and also the value 
of the peak in the circular velocity profile. Both halo finders 
give consistent values for the the final sizes of the subhalo af- 
ter stripping. For the two subhaloes that do not pass through 
the centre of the halo, the amount of stripping is noticeably 
less. The subhalo loses around 35 per cent and 50 per cent of 
its mass for closest approaches of 0.5r v i r and 0.2r v i r respec- 
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Figure 6. The position of the peak of the circular velocity profile 
in relation to the concentrat ion of a halo. Typical halo concentra- 
tions from lNeto et al. ] ( 2007) and radial densities are also labelled. 



tively compared with over 75 per cent if it passes through 
the centre. 

Comparing the halo finders as the subhalo passes 
through the central region of the halo, both show a char- 
acteristic dip in the number of particles recovered. It is 
also noticeable that as the subhalo leaves the centre of the 
halo, AHF always finds a larger subhalo than SUBFIND. This 
is also shown very clearly in Fig. U where in the region 
< r SC p/r v ir < 1 AHF gives much higher recovery of par- 
ticles compared with SUBFIND which has flattened off. The 
cause of this difference can be seen in the lower left panel of 
Fig. [5] by considering the maximum circular velocity. After 
the subhalo has passed through the centre of the halo, the 
maximum circular velocity recovered by AHF spikes meaning 
that background halo particles are being included in the sub- 
halo. There is no such spike in the SUBFIND value (lower right 
panel). This shows that the unbinding of particles is more 
efficient in SUBFIND than AHF. This discrepancy is caused by 
AHF assuming spherical symmetry for the unbinding when 
the subhalo becomes elongated in the centre of halo and is 
no longer a spherical shape. 

For the subhalo with the closest approach of 0.5r v j r , AHF 
shows a smooth transition in the size of the subhalo, while 
SUBFIND shows the size to decrease and then increase again. 
During this transition the subhalo always has a finite size 
as the subhalo does not pass close enough to the halo centre 
to completely vanish. The decrease and increase in the size 
of the subhalo demonstrates that it is being truncated at a 
radius smaller than its actual size. As the sad dle point in the 
densi ty profile corresponds to the tidal radius (jTormen et al.l 
1998), this in turn shows that a subhalo not passing through 
the centre of a halo will not be completely stripped down 
to its tidal radius. This is perhaps not that surprising as 
the subhalo has not spent a long enough time in the halo to 
undergo the full effects of tidal stripping. 



1.6 




0.8 I , , , 

10 100 1000 10000 

N vlr 

Figure 7. The recovered maximum circular velocity compared 
with number of particles used to generate a M V \ T = 10 12 Mq 
and c = 12 halo. Error bars represent 1 standard deviation and 
are distributed symmetrically in log space. For the average to 
be within 2.5 per cent of the maximum value, in excess of 500 
particles are required. 

The maximum circular velocity is shown in Fig. [5] to be 
a much more stable quantity compared to particle number 
as expected from i )3.1l The strong radial dependence of sub- 
find in particle number is not present in maximum circular 
velocity. While this is an advantage in recovering properties 
of the subhalo, Fig. [5] also shows how this quantity can be 
misleading when considering stripping. For the case where 
the subhalo passes within 0.5r v i r , the subhalo was stripped 
of around 35 per cent of its mass, but the maximum circular 
velocity changes by less than 5 per cent. This is caused by the 
maximum circular velocity being located at a radius much 
closer to the centre of the subhalo and so is less affected by 
stripping which occurs primarily in the outer regions. 



5 CIRCULAR VELOCITY 

As seen in the previous sections, the peak in the circular ve- 
locity profile of a subhalo is a more stable quantity to recover 
than the total subhalo mass. The origin of this stability is 
related to the fact that the radius at which the maximum cir- 
cular velocity is reached is located much closer to the centre 
of the halo and so is unaffected by truncation. Fig. [6] shows 
how the position of peak changes with the concentration of 
a halo. For a NFW halo this can be obtained numerically to 
give, 

' vmax 2.16 
rVir c 

The values determined by equation (I14|l are based on an 
ideal NFW halo, but for low resolution haloes there will be 
deviations from this curve. For the subhalo used in this work 
(c = 12) rvmax = 0.18r V ir which corresponds to roughly rsooo 
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(the radius at which the enclosed density is 5,000 times the 
critical density, p C rit)- Stripping occurs in the outer regions 
of the subhalo and so for it to affect this radius a large 
amount of material needs to be lost, consistent with Fig. [5] 
One of the main issues with using the maximum cir- 
cular velocity of a halo is how its measurement depends 
upon resolution. To investigate this, we generated a halo 
with Mvir = 1O 12 M0 and c = 12 in isolation using a differ- 
ent number of particles within the virial radius each time. 
For each number of particles within the virial radius, we 
constructed 1,000 realisations in order to constrain the vari- 
ation. Fig. shows how the recovered maximum circular ve- 
locity varied with the total particle number. For the sparsely 
populated realisations the average maximum circular veloc- 
ity was higher than the analytic value. As more particles 
were used, the two values converged. For the average value 
to be within 2.5 per cent of the analytic value, in excess 
of 500 particles were required in the halo. The variation of 
the maximum circular velocity between different realisations 
of the same total virial particle number is strong for the 
sparsely populated haloes. At all points the curve is within 
f standard deviation of the analytic value, but the variation 
is clear where for fO particles the standard deviation is 0.56 
compared with 0.002 for 10,000. To obtain an accurate value 
for the maximum circular velocity of a recovered subhalo, its 
resolution has to be taken into account. 



6 SUMMARY AND CONCLUSIONS 

Halo finders are an important tool for the analysis of cosmo- 
logical simulations. They are pivotal in the construction of 
merger trees, which underpin galaxy formation modelling, 
and their results allow us to characterise, for example, the 
abundance and spatial distribution of both dark matter 
haloes and subhaloes. There are as many techniques for iden- 
tifying haloes and subhaloes in cosmological simulations as 
there are halo finders and so it is interesting to ask whether 
or not (sub)halo properties recovered by different halo find- 
ers are consistent. 

In this paper we have compared and contrasted the re- 
sults of two halo finders, SUBFIND and AHF, that use funda- 
mentally different approaches to identifying subhaloes. We 
have taken a simple test problem, the identification of a 
NFW subhalo embedded in a more massive NFW halo, and 
compared the performance of SUBFIND and AHF in recovering 
the mass of the subhalo at different radii within its host. As 
shown using SUBFIND, halo finders that identify subhaloes as 
overdensities will have a strong dependence on the local den- 
sity. This is demonstrated in the strong radial dependence 
in the fraction of a model subhalo SUBFIND recovers. As the 
subhalo gets closer to the centre of the halo, the background 
density from the halo is rising. With a higher background 
density and the same density for the subhalo, the overden- 
sity will be less leading to a smaller subhalo being recovered. 
By the time the subhalo is in the centre of the halo, which 
corresponds to the densest point, the overdensity becomes 
negligible leading to no saddle point and the subhalo is no 
longer detected. While the size of the overdensity recovered 
roughly corresponds to the tidal radius of the subhalo, it 
has been shown that not all subhaloes are stripped down 
to this size when they pass through a halo. The authors of 



SUBF IND are aware of these issues (see §4.1 of ISpringel et al.l 
2008) and post-process, but where this effect is not taken 
into account it could have profound consequences on sub- 
structure studies. 

The radial dependence of locating subhaloes as overden- 
sities will have a large effect on measures of tidal stripping. 
As a subhalo plunges into a halo, the halo finder will reduce 
the size of the subhalo due to the increase in density. If this is 
not considered, then it will appear the subhalo is undergoing 
a larger amount of stripping as it falls through the halo than 
it actually underwent. Stripping will be further complicated 
by the fact it occurs in the outer region of the subhalo, an 
area that is not included in the truncated subhalo that is 
recovered. This can lead to confusion when comparing the 
recovery of AHF and SUBFIND. AHF indicates that most of 
the stripping occurs as the subhalo passes through the cen- 
tre of the halo and not during the infall, but AHF has been 
shown to have inefficient unbinding causing it to retain a 
larger fraction of particles. Meanwhile SUBFIND indicates a 
more gradual process, but the effects of truncation will cause 
the recovered subhaloes to always be lower estimates of the 
size. Further studies will need to be made to determine how 
dramatic the effect of stripping is on an infalling subhalo. 

The radial dependence in recovery will also have im- 
portant implications for the subhalo mass distribution. Two 
subhaloes that have identical mass can be recovered with 
different sizes based on position. This will lead to large sub- 
haloes being recovered as smaller ones, in turn, leading to 
subhalo mass distributions biased towards the low mass end. 
Whilst most subhaloes that reside in the inner region of the 
halo will have undergone a large amount of stripping and 
will be smaller anyway, the effect of truncation still needs 
to be considered alongside the underlying physics. These is- 
sues highlight that the recovered mass identified using the 
overdensity method is not a good property to consider when 
studying subhaloes. This is true even as far out as the virial 
radius of the halo, where the mass can be underestimated 
by around 25 per cent. 

A more stable quantity to consider is the peak in the 
circular velocity profile. This is located much closer to the 
centre of the subhalo and so will be less affected by trunca- 
tion and the particular choice of the definition for an entire 
subhalo. Both AHF and SUBFIND recover consistent values 
for the maximum circular velocity at all radii within the 
halo, except at the very centre of the halo where no par- 
ticles are recovered. This makes the circular velocity peak 
a useful quantity to track subhaloes and gives a good indi- 
cation of initial mass. However, when considering stripping, 
the circular velocity peak is no longer useful. Being located 
so close to the centre of the subhalo, a substantial amount 
of the outer layers can be stripped before the peak in the 
circular velocity is affected. 

Two methods of improving the accuracy of subhalo re- 
covery would be halo tracking and phase space. Halo track- 
ing involves identifying the subhalo before it falls into the 
halo so all the particles that were originally part of the struc- 
ture are followed and at each time step they can be tested to 
see if they are still part of the substructure. The disadvan- 
tage of this technique is that it requires multiple snapshots 
to identify the subhalo, not a problem for the second method 
of phase space. Phase space takes into account not only the 
spacial position of the subhalo particles, but also links par- 
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tides based on a common velocity. By considering haloes 
in phase space density, any subhaloes that are present will 
stand out as overdensities. These can then be isolated. For 
subhaloes in the centre of the halo, the difference in the 
bulk velocity of the particles would cause them to be sep- 
arated in phase space. The only remaining problem would 
be if a subhalo was at rest in the centre of the halo. These 
structures could not be separated in phase space, but it is 
arguable whether such a structure would be a dynamically 
independent entity. 
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